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For  a  mathematical  model  based  on  the  result  of  physical  measurements,  it  becomes  possible  to 
determine  their  influence  on  the  final  solution  and  its  accuracy.  However,  in  classical  approaches,  the 
influence  of  different  model  simplifications  on  the  reliability  of  the  obtained  results  are  usually  not 
comprehensively  discussed.  This  paper  presents  a  novel  approach  to  the  study  of  methane/steam 
reforming  kinetics  based  on  an  advanced  methodology  called  the  Orthogonal  Least  Squares  method.  The 
kinetics  of  the  reforming  process  published  earlier  are  divergent  among  themselves.  To  obtain  the  most 
probable  values  of  kinetic  parameters  and  enable  direct  and  objective  model  verification,  an  appropriate 
calculation  procedure  needs  to  be  proposed.  The  applied  Generalized  Least  Squares  (GLS)  method  in¬ 
cludes  all  the  experimental  results  into  the  mathematical  model  which  becomes  internally  contradicted, 
as  the  number  of  equations  is  greater  than  number  of  unknown  variables.  The  GLS  method  is  adopted  to 
select  the  most  probable  values  of  results  and  simultaneously  determine  the  uncertainty  coupled  with  all 
the  variables  in  the  system.  In  this  paper,  the  evaluation  of  the  reaction  rate  after  the  pre-determination 
of  the  reaction  rate,  which  was  made  by  preliminary  calculation  based  on  the  obtained  experimental 
results  over  a  Nickel/Yttria-stabilized  Zirconia  catalyst,  was  performed. 

©  2014  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

The  methane/steam  reforming  process,  often  used  for  the  syn¬ 
thesis  gas  (syngas)  production  from  hydrocarbon  source,  is  one  of 
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the  fuel  conversion  methods,  which  are  especially  important  in  the 
context  of  high  temperature  fuel  cells  such  as  Solid  Oxide  Fuel  Cells 
(SOFCs).  This  technology  is  beneficial  for  power  generation  in 
SOFCs,  considering  their  fuel  flexibility  [1  ]. 

The  high  operating  temperature  of  the  cell  eliminates  the  ne¬ 
cessity  of  the  external  heat  source  for  the  strongly  endothermic 
reforming  reaction.  Regarding  this  advantage,  coupling  of  the 
reforming  process  with  SOFCs  is  a  beneficial  solution  for  hydrogen 
production  in  fuel  cell  technology.  To  properly  address  and  imple¬ 
ment  this  technique  to  the  technological  development  of  SOFCs,  the 
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implementation  of  an  internal  reforming  system,  modelling  and 
numerical  analyses  of  the  phenomena  occurring  inside  the  SOFC 
systems  are  important  subjects.  The  reason  relies  on  the  fact  that 
the  fuel  reforming  process,  in  particular  the  designing  and  opti¬ 
mizing  of  an  indirect  internal  reformer  and  a  cell  with  direct 
reforming  for  an  SOFC  application,  involve  problems  connected 
with  heat  transfer  phenomena  in  porous  media.  To  conduct  nu¬ 
merical  modelling,  proper  kinetics  of  the  reaction  is  also  needed  to 
improve  the  reliability  of  numerical  modelling.  Thus,  the  treatment 
of  reaction  kinetics  and  physical  phenomena  in  a  system  are  very 
important  subjects  to  be  accomplished. 

The  kinetics  of  the  reforming  process  can  be  described  by  one  of 
three  concepts:  i)  General  Langmuir-Hinshelwood  kinetics,  ii)  first 
order  reaction  with  respect  to  methane  and  iii)  power  law  ex¬ 
pressions  derived  from  data  fitting  [2—5],  One  of  commonly  used 
mathematical  expressions  of  the  methane/steam  reforming  reac¬ 
tion  is  determined  as  the  reaction  rate:  rCH4  =  k(PcH4)a(PH2o),J.  by 
taking  the  advantage  of  the  power  law  expression.  The  reaction  rate 
of  the  reforming  reaction  is  very  sensitive  on  the  catalytic  activity, 
reaction  temperature  and  the  partial  pressure  of  the  reactants. 
Thus,  the  investigation  relies  on  the  material  properties  and 
transport  phenomena.  To  find  out  the  proper  reaction  kinetics,  the 
preparation  of  the  sample  and  experimental  conditions  are  very 
important  concerns  to  clarify  the  reaction  kinetics. 

There  are  many  studies  focusing  on  the  kinetics  of  methane/ 
steam  reforming  [2—16],  The  experimental  investigation  has  been 
also  done  over  the  Nickel/Yttria-Stabilized  Zirconia  (Ni/YSZ),  which 
is  a  popularly  adopted  SOFC  anode  material  [2,3,6-13,15,16],  As 
seen  in  Refs.  [2,3,6—13,15,16],  the  derived  empirical  parameters 
describing  the  process  reported  in  them  are  divergent  among 
themselves.  For  instance,  Table  1  shows  one  of  empirical  parame¬ 
ters,  the  so-called  reaction  order  a  and  b  found  in  the  literatures 
[2,7,9-13,15,16],  Although,  the  empirical  parameters  a  and  b  should 
be  independent  from  the  specific  property  of  the  investigated 
sample,  as  far  as  the  same  type  of  catalyst  is  used  [9],  Moreover, 
uncertainties  of  the  obtained  result,  which  may  have  influenced  the 
derived  reaction  kinetics,  were  not  substantially  discussed  in  the 
published  investigations.  Thus,  this  paper  presents  a  new  approach 
to  the  analysis  of  the  reforming  process,  introducing  a  numerical 
procedure  incorporating  the  Generalized  Least  Squares  (GLS) 
method. 

The  conventional  methods  of  modelling  multi-physical  phe¬ 
nomena  involving  (heat  and  mass)  transport  problems  usually 
describe  analysed  systems  with  a  set  of  equations,  which  are 
uniquely  defined  in  the  mathematical  sense.  In  general,  this  type  of 
theoretical  approach  yields  unique  solutions  from  the  peculiarly 
specified  constraints  set  of  equation  for  a  system.  Besides,  they 
assume  that  each  variable  or  model  parameter  is  precisely  defined 
and  their  uncertainties  are  negligible  (Fig.  1(a)).  However,  in  the 
case  of  the  methane/steam  reforming  process,  in  which  the  model 
parameters  include  the  properties  from  a  general  thermodynamic 


parameter,  reaction  orders  a  and  b,  found  in 


Reference  a 

Brus  [2]  0.98 

Ahmed  and  Foger  [7]  0.85 

Iwai  et  al.  [9]  0.82 

King  et  al.  [10]  1 

Leeetal.  [11]  1 

Odegard  etal.  [12]  1.2 

Timmermann  et  al.  [13]  1.19 

Yakabeetal.  [15]  13 

Mogensen  [16]  0.7 


-0.09 

-035 

0.14 

0 

-1.25 

0 

0 
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table,  the  assumption  of  accuracy  of  all  of  those  factors  can  lead  to 
the  simplification  of  the  model.  The  model  simplification  can  also 
be  inducted  by  an  applied  mathematical  approach  both  in  the 
phase  of  definition  and  during  the  model  solving  stage.  Moreover, 
one  of  the  most  significant  factors  in  identifying  the  inaccuracies  of 
the  model  solutions  is  connected  with  the  uncertainties  of  the 
directly  measured  variables  (observations)  incorporated  in  the 
problem  description  (Fig.  1(b)).  The  model  defined  in  this  manner 
results  in  the  unique  description  of  the  phenomena,  although  the 
numerical  results  can  differ  from  the  actual  state  of  the  system 
behaviour  and  the  degree  of  accuracy  of  the  model  that  is  usually 
difficult  to  estimate.  The  present  study  analyses  the  concept  of  most 
probable  value  of  a  theoretical  solution  and  a  measure  of  its  inac¬ 
curacy.  The  proposed  methodology  allows  for  including  additional 
data  directly  in  the  mathematical  model  of  the  process  as  the  new 
(supplementary)  variables.  A  model,  which  includes  supplemen¬ 
tary  data,  contains  more  equations  than  unknowns,  and  because 
the  uncertainties  of  the  directly  measured  variables  (observations) 
can  be  internally  contradicted.  Choosing  different  subsets  of 
constraint  equations  allows  for  the  calculation  of  the  unknowns  in 
many,  theoretically  correct  ways  and  a  problem,  which  arises,  is 
how  to  determine  the  most  probable  result.  The  GLS  procedure  is 
proposed  as  a  method  of  correcting  the  measured  data,  securing  its 
higher  accuracy  and  obtaining  the  most  probable  value  and  its 
uncertainty  for  the  parameters  to  be  determined  by  numerical 
computation  (Fig.  1(c)).  Additionally,  it  can  be  noted  that  the  GLS 
method  provides  the  objective  criteria  for  the  formal  evaluation 
and  falsification  of  different  mathematical  models  of  investigated 
phenomena. 

The  present  paper  introduces  a  novel  approach  to  derive  reac¬ 
tion  kinetics  described  by  a  power  law  type  of  equation.  The  effect 
of  uncertainties  on  the  derived  reaction  rate  is  primarily  discussed 
by  adopting  the  GLS  method.  This  approach  leads  to  improving  the 
reliability  of  the  derived  reaction  kinetics  of  the  methane/steam 
reforming  process.  The  challenge  to  conclude  the  applicability  of 
the  GLS  method  in  the  analysis  of  the  chemical  reaction  process  is 
also  highlighted  in  this  paper. 

2.  Computational  method 

The  classical  approach  to  mathematical  modelling  assumes  that 
real  life  problems  can  be  described  and  defined  by  a  finite  set  of 
mathematical  equations.  In  general,  the  problem  with  N  unknowns 
and  K  measured  values  can  be  solved  unambiguously  when  there 
are  N  equations  in  the  mathematical  model.  However,  in  many 
theoretical  approaches,  the  uncertainties  created  by  model  sim¬ 
plifications  and  experimental  errors  cause  the  unknown  uncer¬ 
tainty  of  the  obtained  solution.  On  the  other  hand,  when  the 
number  of  equations  exceeds  the  number  of  unknowns  (a  set  of 
overdetermined  equations)  and  the  measured  values  are  charac¬ 
terized  by  experimental  errors,  the  model  equation  set  becomes 
internally  contradicted.  The  Least  Squares  (LS)  method  provides  a 
strategy  for  finding  the  approximate  and  most  probable  solution  of 
an  overdetermined  system  of  equations.  The  first  comprehensive 
formulation  of  the  least  squares  (LS)  method,  was  presented  by 
Legendre’s  in  1805  [17],  Initially,  the  LS  method  was  used  mainly  in 
geodetic  calculations  and  then  Sweneker,  who  investigated  the 
balances  of  substances  in  the  branched  pipelines,  presented  its  first 
application  to  the  energy  science  problem  [18],  The  application  to 
the  energy  science  and  the  theoretical  aspects  of  the  LS  method 
were  developed  [19-26],  The  research  involved  problems  con¬ 
nected  with  theory  of  coordination  of  material  and  energy  balance 
[19—21],  heat  transfer  by  convection  and  radiation  and  the  pre¬ 
diction  of  temperature  distribution  [22,23],  Other  concerned  issues 
were  investigations  related  to  carbon  dioxide  emissivity  [24], 
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Fig.  1.  Approaches  to  mathematical  modelling:  (a)  classical  problem  definition,  (b)  errors  of  mathematical  models  and  (c)  models  with  supplementary  data. 


solidification  of  binary  alloys  [25],  and  modelling  of  the  autocata- 
lytic  reactions  [26], 

2.1.  Least  squares  method 

Mathematical  models  describing  real  process  are  usually 
defined  by  nonlinear  equation  sets  consisting  of  algebraic  functions 
and  differential  equations.  First,  assume  that  the  considered  prob¬ 
lem  is  described  by  a  set  of  J  algebraic  equations: 

fuA  (x\\ 


Moreover,  it  is  assumed  that  experimental  errors  Sk  of  mea¬ 
surements  satisfies  the  normal  distribution  with  mean  value  equals 
zero  and  variances  a'l :  sk~N(0,ofy. 

According  to  assumptions,  fj  are  the  unrestricted  algebraic 
functions.  However,  the  least  squares  method  requires  the  linear 
form  of  the  constraint  equations.  If  the  fj  functions  are  differen¬ 
tiable,  they  can  be  presented  in  the  form  of  a  Taylor  series  in  the 
neighbourhood  of  the  point  P(u,x).  The  linearization  of  a  function  fj 
is  the  first  order  term  of  its  Taylor  expansion  around  the  point  of 
interest.  After  linearization,  the  constraint  equations  can  be  written 
in  the  matrix  form: 


=  0,  0  =  1,2,  ,J),  u  = 


[4] 


;*  y  where  A  =  dfi/duk  is  J  x  K  Jacobi’s  matrix  and  B  =  d/j/dx„  is  J  x  N 
,  Jacobi’s  matrix.  The  LS  method  assumes  the  minimization  of  the 
^  '  following  function: 


where  u  is  K-elements  vector  of  the  quantities  to  be  measured  and 
x  is  N-elements  vector  of  the  quantities  to  be  solved  by  numerical 
computation.  Although  for  measured  values  u  and  approximations 
of  unknowns  x  the  equation  set  is  not  precisely  fulfilled: 


(2) 


m  **  vJCAv  -  (5) 

where  Cs  is  the  covariance  matrix  with  the  measurement  variances 
al  on  the  diagonal  and  0  on  the  other  positions  (due  to  the  inde¬ 
pendence  of  the  variables).  The  method  of  Lagrange  multipliers  is 
applied  in  order  to  achieve  the  extreme  and  find  corrections: 

y  =  £HbtF-1w  (g) 

v  |=  CsAJF-'(w-By) 


where  Wj  designates  the  residual  of  the  j-th  equation.  Real  values  u 
and  x  are  sums  of  measurements  and  approximations  of  un¬ 
knowns  (u  and  x)  and  corrections  (v  and  y):  u*  =  u  +  v  and 
x*  =  x  +  y.  Therefore  it  can  be  written: 


where  F  =  ACsAT  and  G  =  BTF  'B. 

2.2.  Generalized  Least  Squares  method 


fj(u  +  v,x+y)  =  0,  (jj  —  1.2. •••,/), 


(3) 


The  Generalized  Least  Squares  (GLS)  method  (Unified  Least 
Squares  Method)  was  introduced  by  Mikhail  and  Ackermann  [27], 
The  classical  definition  of  the  least  squares  method  assumes  that 
the  rank  of  matrix  A  is  equal  to  the  amount  of  model  equations  J.  In 
many  problems  this  assumption  is  not  fulfilled  because  the  number 
of  equations  fj  is  different  than  the  amount  of  measured  values  Uk. 
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Thus,  Eq.  (6)  proposed  in  the  classical  method  cannot  be  used 
because  the  matrix  F-1  does  not  exist. 

The  most  important  assumption  of  the  GLS  method  is  that  all 
model  variables  (measured  values  u  and  unknowns  x)  are  obser¬ 
vations  and  treated  numerically  in  the  same  manner.  The  only 
difference  between  them  is  the  assumption  that  errors  of  un¬ 
knowns  are  much  bigger  than  the  errors  of  measured  values 
(sk»s„). 

After  the  linearization  of  the  model  equations  the  constraints 
can  be  written: 


(7) 


where  Ab  =  [AB]and  Vb  —  [vy]Tand  then  the  covariance  matrix  has 
the  following  form: 


CB 


0  1 

CSXJ 


(8) 


where  Cs  and  Csx  are  the  covariance  matrices  with  diagonals  con¬ 
sisting  of  variances  of  measurements  <j\  and  unknowns  <r„ 
respectively.  To  solve  the  GLS  problem  and  minimize  the  function 
$(v,y),  which  has  the  following  form: 

*(v,3 o  =  E  d)2  +  E  (S)2^min  (9) 


the  Lagrange  multiplier  method  is  used  and  it  gives  the  solution: 

VB  =  CbA^Fb'Wb,  where  FB  =  AbCba£  (10) 


2.3.  Errors  of  solution 


To  yield  new  error  bounds  for  corrected  measured  data  (u,  Cu) 
and  the  theoretically  solved  variables  (x,  Cx),  the  law  of  error 
propagation  can  be  applied.  The  new  covariance  matrices  for  the 
classical  LS  method  can  be  written: 


Cu  =  Cs-SACs 
Cx  =  (BtT_1B)  1 


(11) 


where  S  =  CsA'f'/IBG  'b'F  ’)  and  I  is  the  unit  matrix. 

On  the  other  hand,  in  the  case  of  the  generalized  approach,  the 
calculated  covariance  matrix  of  improved  observations  has  the 
following  form: 


Cvb  =  Cb-  CbAtbFb'AbCb  =  Cc^] 


(12) 


2.4.  Interpretation  of  the  covariance  matrix 


The  covariance  matrix  Cvb  defined  by  Kolenda  et  al.  [25]  rep¬ 
resents  the  hyperellipsoid  of  the  (K  +  N)-dimensional  normal  dis¬ 
tribution.  The  covariance  ellipsoid  is  the  hypersurface  in  the 
(K  +  N)-dimensional  space  and  is  always  inside  the  (K  +  N)- 
dimensional  cuboid  defined  by  the  central  point  P'(u  ,x  )  and  the 
square  roots  of  the  diagonal  elements  of  Cvb  which  represent  the 
lengths  of  the  cuboid  edges.  Covariance  hyperellipsoid  touches 
cubicoid  in  2(I<  +  N)  points  and  its  volume  can  be  calculated  from 
the  formula: 


VK+N  = 


\  iK  i  N  ,  2Wd<'N.2  |CvB 
fQ(I<  +  N)  +  l) 


where  T  is  the  Euler’s  gamma  function. 


(13) 


Generally,  it  can  be  noted  that  the  lines  of  constant  probability 
are  ellipsoids  similar  to  the  covariance  ellipsoid.  Those  ellipsoids 
are  situated  inside  (outside)  of  it  for  larger  (smaller)  probability.  For 
the  largest  probability  the  ellipsoid  collapses  to  the  central  point 
[28],  Thus,  it  can  be  pointed  out  that  conditions  for  achieving  the 
global  minimum  can  be  presented  equivalently  as  the  minimization 
of  the  hyperellipsoid’s  volume  and  the  minimization  of  the 
covariance  matrix’s  determinant: 

min(l/K ,  N)  ->min(|CvB|)or  min  ^  E  +  E  an 

(14) 

The  presented  formula  can  be  used  for  the  falsification  of 
different  mathematical  models  describing  the  examined  physical 
phenomena. 


3.  Experimental  studies  on  methane/steam  reforming 

The  analysis  of  the  reforming  process  was  coupled  with  the 
experimentation  of  methane/steam  reforming  reaction  kinetics 
derived  over  a  Nickel/Yttria-Stabilized  Zirconia  catalyst  material 
(NiO/YSZ  60:40  vol%,  AGC  SEIMI  CHEMICAL  CO.  LTD  [29],  The  Ni / 
YSZ  cermet  is  chosen  for  the  anode  of  SOFCs,  because  this  material 
satisfies  the  following  criteria  of  requirements:  inexpensive  fabri¬ 
cation  cost,  high  electrical  conductivity,  high  electrocatalytic  ac¬ 
tivity,  stability  under  reducing  environment,  thermal  expansion 
and  chemical  compatibility  with  other  materials  composing  the 
SOFC  cell  layers  and  ensuring  sufficient  Three  Phase  Boundary 
(TPB)  for  the  electrochemical  process  [30—33].  Moreover,  a  lot  of 
studies  on  the  reaction  kinetics  of  the  methane/steam  reforming 
reaction  over  the  Ni/YSZ  catalyst  were  published  earlier  and  huge 
inconstancies  were  found  in  them  as  preceded  in  the  introduction. 
Thus,  the  main  goal  of  the  present  work  is  to  evaluate  the  effect  of 
acceptable  inaccuracies  (caused  by  the  specificity  of  measure¬ 
ments)  on  the  reforming  performance  in  the  form  of  the  rate 
equation,  which  is  dependent  on  the  empirical  parameters.  To 
satisfy  the  objectives,  the  combined  experimental-numerical 
methodology  proposed  by  Odegard  [12]  and  modified  by  Brus 
et  al.  [8]  was  applied. 

Fig.  2(a)  shows  the  schema  of  the  experimental  setup.  Three 
main  parts  of  the  system  can  be  noted:  the  inlet  part,  where  the 
substrate  gases  are  prepared,  the  reaction  part  with  the  reformer 
tube  placed  in  the  furnace  and  the  outlet  part  where  the  product  of 
the  reforming  reaction  is  analysed. 

The  experiment  was  performed  with  high  purity  methane  (CH4) 
supplied  from  a  gas  cylinder  via  the  mass  flow  controller  and  water 
(H2O)  fed  with  the  pump.  Additionally,  nitrogen  (N2)  was  also 
supplied  to  the  reactor  to  maintain  the  reforming  conversion  rate  at 
a  low  level  and  to  enable  to  derive  the  correct  reforming  kinetics 
equation,  since  the  reaction  has  to  occur  in  the  entire  volume  of 
catalyst.  N2  does  not  have  a  direct  influence  on  the  occurring  re¬ 
actions  but  changes  the  partial  pressures  of  the  components  and 
finally  decreases  the  methane  conversion  rate.  Thus,  the  gas 
mixture  of  CH4,  N2  and  H20  was  supplied  to  the  reaction  zone.  CH4 
was  regulated  200  [kPa]  and  supplied  with  the  mass  flow 
controller,  N2  was  and  also  treated  in  the  same  way.  H2O  was 
pumped  and  supplied  to  the  evaporator,  which  is  used  at  the  same 
time  as  the  preheater  for  heating  the  gas  mixture. 

The  main  part  of  the  system  is  the  stainless  steel  reformer  of 
which  the  schema  is  presented  in  Fig.  2(b).  In  the  experimental 
investigation  the  Nickel/Yttria-Stabilized  Zirconia  catalyst  in  the 
form  of  fine  quality  powder  of  NiO/YSZ,  of  which  the  composition  is 
60:40  in  volume  ratio,  was  used.  The  catalyst  powder  used  in  the 
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Fig.  2.  Schema  of  (a)  experimental  setup  and  (b)  reactor. 


experiment  is  spherical  in  shape  and  has  a  0.85  [pm]  diameter  with 
a  specific  surface  area  5.2  [m2]  per  1  [g]  of  catalyst  [29],  The  fine 
powder  catalyst  was  used  in  order  to  determine  the  kinetics  of  the 
catalysed  methane/steam  reforming  process  without  the  influence 
of  mass  transport  process  [16],  It  is  crucial  to  avoid  the  mass 
transfer  limitations  because  they  can  lower  the  activity  of  the 
catalyst,  especially  in  the  case  of  exponentially  fast  reactions  [Il¬ 
ls, 16],  The  catalyst,  of  which  the  amount  used  in  the  experiment 
was  3  [g],  was  dispersed  to  be  a  homogeneous  distribution  over  the 
quartz  wool-made  bed  located  in  the  middle  of  the  reactor.  Then, 
the  catalyst  layer  was  covered  with  the  quartz  wool  and  then  the 
half  upper  of  the  reactor  was  fulfilled  with  AI2O3  balls.  The  catalyst 
is  preceded  by  AI2O3  balls  used  to  prevent  a  cooling  effect  of  the 
entering  fluid  and  avoid  a  large  temperature  gradient  in  the 
reformer.  The  gas  mixture  reaches  the  reaction  zone  with  the 
catalyst  after  preheating  by  the  electric  furnace  to  the  reaction 
temperature.  A  reformer  tube  is  placed  in  the  electrical  furnace 
with  a  maximum  temperature  of  800  [°C[.  The  preheater  and  after¬ 
heater  allow  to  reach  the  maximal  temperature  of  400  [°C],  how¬ 
ever,  the  temperature  was  set  at  200  [°C]  for  both  of  them  during 
the  experimental  procedure.  The  temperature  of  the  reactor  and 
the  temperature  around  it  are  controlled  with  the  thermocouples 
located  before  and  after  the  furnace  as  well  as  by  the  thermocou¬ 
ples  vicinity  of  the  reformer  (on  the  surface  and  inside  reaction 
zone). 

After  the  reforming  reaction  the  exhaust  gas  was  analysed  by 
the  gas  chromatograph  (GC390B  by  GL  Sciences).  The  gas  supplied 
to  the  gas  chromatograph  was  dried  by  cooling  down  the  gas 
mixture  to  2  [°C]  for  the  preliminary  of  the  chromatography.  Two 
methods  of  gas  detection  were  used:  TCD  (Thermal  Conductivity 
Detector)  to  the  methane  m“H)  and  hydrogen  and  FID  (Flame 
Ionization  Detector)  with  a  methanizer  for  the  carbon  monoxide 
W.QQ  and  carbon  dioxide  m°0  detection,  where  m  stands  for  the 
molar  fraction  of  the  chemical  species  and  the  superscript  0  de¬ 
notes  measurements  at  the  outlet. 

Before  conducting  measurements,  the  catalyst  was  treated  at 
the  evaluated  temperature  of  800  [oC]  with  a  mixture  of  150 
[ml  min-1]  nitrogen  and  100  [ml  min-1]  hydrogen  due  to  the 
reduction  process  of  NiO.  NiO/YSZ  had  been  kept  in  the  reduction 
process  for  6  [h].  The  measurements  were  carried  out  with  different 
conditions  based  on  the  various  temperatures  (in  the  range  of  550— 
750  [°C]),  Steam-to-Carbon  ratio  SC  (in  the  range  of  2.5—6)  and 
Nitrogen-to-Carbon  ratio  NC  (in  the  range  of  0.5—6).  The  details  of 
the  experimental  conditions  are  shown  in  Table  2.  The  other 


perspective  of  supplying  N2  was  for  levelling  the  residential  time  of 
the  mixture  gas  flow  at  constant  among  the  various  flow  rate 
conditions. 


4.  Analysis 

4.1.  Preliminary  analysis 


Methane/steam  reforming  is  the  conventional 
hydrogen  production  and  reactions  governing  this 
following  two  reactions  [12  !. 

Methane/steam  reforming  reaction: 

method  of 
process  are 

CH4+H2O— -3H2  +CO 

(15) 

Shift  reaction: 

CO  +  h2o~h2+co2 

(16) 

The  reaction  rate  of  the  methane/steam  reforming  reaction  can 
be  described  by  the  following  equations: 

r  =  fc(PcH4)a(PH2o)b 
k  =  AexpCRj) 

(17) 

where  r  is  the  reaction  rate  [mol  s-1  g-1  ],  k  is  the  reaction  constant, 
Pch4  and  Ph2o  are  the  partial  pressures  of  CH4  and  H2O  [Pa],  a  and  b 
are  the  reaction  order  coefficients  [-],  A  is  the  pre-exponential 
factor  [mol  g-1  s-1  Pa-(a  +  b)],  E  is  the  activation  energy  0  mol-1], 
R  is  the  universal  gas  constant  [J  mol-1  K-1]  and  T  is  the 


Table  2 

Experimental  conditions  for  different  values  of  SC  and  NC  ratios. 


SC  [-]  NC  [-]  Volume 

[ml  min-1] 
CH4 

3  3  50 

2.5  6  37 

5  4  35 

4  2  50 

3  1  70 

6  0.5  47 


Volume  Volume 

flow  rate  flow  rate 

[ml  min"1]  [ml  min"1] 

N2  h2o 

150  0.11 

221  0.07 

140  0.13 

100  0.15 

70  0.15 

23  0.21 


1.43  x 
1.43  x 
1.43  x 
1.43  x 
1.43  x 
1.43  x 


10~2 

10"2 

10~2 
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temperature  [K],  On  the  other  hand,  the  shift  reaction  is  a  fast 
process  and  reaches  equilibrium  instantaneously,  and  therefore  it  is 
assumed  that  its  conversion  rate  can  be  calculated  from  the  equi¬ 
librium  equation: 


Ph2PC02 

PC0Ph20 


(18) 


where  pHl,  pC02,  pc o  and  pH2o  are  the  partial  pressures  of  H2,  C02, 
CO  and  H20  [Pa]  and  AG  is  the  change  in  Gibbs  free  energy 
LI  mol-1]. 

The  reforming  process  is  defined  by  the  values  of  the  empirical 
parameters  describing  the  reaction  order  coefficients  a  and  b, 
activation  energy  E,  and  pre-exponential  factor  A.  To  determine 
their  values  and  the  uncertainty  of  the  obtained  results  the  GLS 
method  has  been  applied.  Two  main  phases  of  the  analysis  can  be 
pointed:  preliminary  calculations  to  estimate  the  initial  values  of 
the  unknown  parameters  in  the  vector  x  and  the  application  of  the 
GLS  method.  The  GLS  algorithm,  which  is  presented  in  Fig.  3, 
started  from  the  determination  of  the  constraint  equations 
describing  the  problem  and  definition  of  the  vector  that  contains 
measured  and  unknown  variables.  The  next  step  was  to  introduce 
the  covariance  matrix  containing  the  uncertainties  of  all  the 
involved  variables.  After  that,  the  values  of  the  unknowns  were 
initialized,  based  on  the  preliminary  calculation.  The  linearization 
of  the  constraint  equations  with  the  determination  of  the  Jacobi’s 
matrix  and  Lagrange  multipliers  method  allowed  for  the  calcula¬ 
tion  of  the  correction  vectors.  In  addition,  new  error  bounds  for 
corrected  measured  data  were  estimated.  The  last  step  was  to  verify 
the  convergence  of  the  obtained  solution  and  stipulate  the  criteria 
of  convergence  or  recalculation  on  the  basis  of  the  unknown  values 
determined  in  the  previous  step. 

The  preliminary  calculations  were  divided  into  two  parts.  The 
first  of  them  was  based  on  the  observation  that  the  value  of  the 
reaction  constant  does  not  depend  on  the  SC  and  NC  ratios  [12],  The 
computation  was  fist  conducted  to  derive  the  reaction  orders,  a  and 
b,  with  the  carefully  designed  experimental  conditions.  Moreover, 
the  SC  and  NC  ratios  characterizing  the  experimental  conditions 
were  widely  and  properly  dispersed  to  minimize  the  influence  of 
the  possible  measurement  uncertainties.  Different  combinations  of 
the  parameters  a  (range  0—1.5)  and  b  (range  —0.5— 0.5)  were  tested 
and  the  calculated  reaction  constants  were  compared  with  the 
experimental  results  from  the  conditions  (Table  2),  at  which  widely 
dispersed  SC  and  NC  ratios  were  applied.  The  range  of  parameters  a 
and  b  was  chosen  in  accordance  with  published  data  [8,10,12],  The 
experimental  conditions  were  set  to  keep  the  total  inlet  flow  rate  at 
a  constant  level.  For  each  computation  conducted  for  different 
experimental  conditions,  the  parameters  a  and  b  with  results 
showing  the  smallest  standard  deviation  from  k  were  selected.  The 
relative  standard  deviation  RSD/<  is  expressed  as  follows: 


t/Ej-t  -  kave)2 

RSDfc  =  v  J  ' - _  (i=l,2,- 


where,  n  and  fcave  stand  for  the  number  of  measurements  and 
average  of  the  reaction  constant  k  at  each  temperature  condition, 
respectively.  In  this  stage,  the  reaction  constant  k  was  derived  by 
the  following  observation:  the  non-equilibrium  reaction  for  the 
plug-flow  reactor  rate  can  be  described  as  a  ratio  between  the 
change  in  the  flow  rate  of  methane  and  the  change  in  the  amount  of 
catalyst  [34],  Mathematically,  the  reaction  rate  of  the  steam 
reforming  reaction  can  be  described  as: 


(20) 


where,  FCH4  =  F‘-(l-x),F‘  is  the  methane  flow  rate  [mol  s  1  ] 
and  one  at  the  inlet  and  of  the  catalyst  layer,  x  is  the  dimensionless 
converted  amount  of  methane  in  the  catalyst  layer,  and  wcat  is  the 
weight  [g]  of  catalyst.  By  combining  Eqs.  (17)  and  (20),  the  reaction 
constant  k  is  given  as  the  following  formula: 


- 1 - J-dx  (21) 

V"W  J  [(PmYiPwf} 

From  the  stoichiometry  of  the  reaction  [15]  and  [16],  the  partial 
pressures  can  be  defined  as: 


Ph-°=  (%f)'P‘T 


1  +  SC  +  JVC  +  2x 
SC-x-v 


(22) 

(23) 


where,  Fan  is  the  total  amount  of  product  [mol  s-1]  at  the  outlet  of 
the  reactor,  P  is  the  total  pressure  [Pa],  Fch4,  Fh2o  are,  respectively, 
the  amount  of  CH4  and  steam  in  [mol  s  ']  at  the  reformer  output 
and  y  is  a  fraction  of  reacted  carbon  monoxide.  Introducing  Eqs. 
(22)  and  (23)  into  Eq.  (21),  the  final  formulation  of  the  reaction 
constant  becomes  as  following: 


fc = T  [  (i+sc+Nc+2*>o+bj  .dx 

V'W  J  |pa+f>(i  -  x)a(SC  -x-  y)hJ 


(24) 


where  a  and  b  are  coefficients  related  to  the  reaction  order.  The 
methane  conversion  rate  x,  which  was  used  to  introduced  to  Eq. 
(24),  is  determined  by  the  outlet  quantity  of  CH4,  CO  and  C02  and 
expressed  as  [35]: 


(25) 


where,  men,’  mco’  mco  's  t^e  mole  fraction  of  CFLt  and  CO  and  C02 
at  the  outlet  of  the  reactor.  The  values  a  =  0.97  [— ]  and  b  =  -0.08  [— ] 
were  characterized  by  the  smallest  summed  standard  deviation  and 
chosen  as  the  best  description  of  the  process  (Fig.  4). 


Fig.  3.  GLS  algorithm  diagram. 
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The  next  stage  was  the  preparation  of  the  Arrhenius  plot  based 
on  the  data  obtained  from  the  experiments  carried  out  in  different 
temperatures  (550-750  [°C]).  The  approximation  line  (Fig.  5): 

lnfe  =  lnA-|i  (26) 

was  used  to  find  the  values  of  the  activation  energy 
E  =  117.219  x  103  [J  mol-1]  and  the  pre-exponential  factor 
A  =  1.554  x  10-3  [mol  g-1  s-1  Pa-0  89].  Thus,  the  first  approximation 
of  the  kinetic  equation  of  the  methane/steam  reforming  reaction 
can  be  written  as: 

Rst  =  Wcat-r  =  Wear  1.554  x  10  3 

/-117xl03\,  ,0.97  /  ,-0.08  (27) 

•  eXP  [ - gf - )  (PcH„)  (Ph2o) 

where,  wcat  stands  for  the  weight  of  a  catalyst  [g]  and  the  subscript 
st  denotes  for  the  steam  reforming.  Fig.  6(a)  shows  the  represen¬ 
tative  examples  of  conversion  rates  obtained  from  the  experiment. 
The  one  of  the  most  important  factors  in  deriving  the  proper 
non-equilibrium  kinetic  equation  are  the  adequate  experimental 
conditions  which  ensure  the  conversion  rate  and  reaction  rate  are 
not  following  the  equilibrium  reaction.  This  can  be  achieved  by 
choosing  the  appropriate  amount  of  catalyst  and  flow  rates  of  the 
reaction  substrates.  Fig.  6(b)  and  (c)  represent  Gas  Hourly  Space 
Velocity  (GHSV)  connecting  to  both  of  those  factors  in  correlation 
with  achieved  conversion  rates.  The  GHSV  is  described  as  the 
function  of  the  mixture  gas  velocity  Vgas  [  mm3  h  1  ]  and  the  volume 
of  the  catalyst  layer  Vcat  [mm3]: 

GHSV  =  (28) 


Fig.  6(b)  is  the  case  where  the  SC  and  NC  ratios  are  fixed  at  3  for 
each  of  them.  In  the  meantime,  Fig.  6(c)  is  the  case  where  the 
temperature  is  fixed  at  700[°C],  It  can  be  noted  that  for  small  GHSV 
values  (small  amount  of  gas  in  high  amount  of  the  catalyst),  the 
achieved  conversion  rate  are  high,  which  means  that  the  reaction  is 
very  close  to  equilibrium  and  is  not  proper  for  delivering  kinetic 
equation.  On  the  other  hand,  for  big  GHSV  values  the  changes  in 
obtained  conversion  rates  are  very  small  and  difficult  to  observe. 
For  kinetic  consideration  the  best  are  moderate  values  of  GHSV, 
which  results  in  moderate  conversion  rates  for  different  observed 
temperatures  and  SC  and  NC  conditions. 


Fig.  4.  Preliminary  calculations  of  the  process  parameters:  analysis  of  the  relative  standard  deviation  (RSD)  of  the  value  of  reaction  constant  at  the  reaction  temperature  of  (a)  700, 
(b)  725  and  (c)  750  [°C]  and  (d)  the  sum  of  the  Relative  Standard  Deviation  (RSDsum)  with  different  values  of  parameters  of  reaction  orders  a  and  b. 
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4.2.  Analysis  with  the  application  of  the  GLS  method 

The  constraint  equations  used  in  the  GLS  algorithm  to  describe 
the  process  of  methane  reforming  can  be  divided  into  three  basic 
groups:  the  kinetic  reforming  equation  (29),  the  linear  Arrhenius 
equation  (30)  and  the  balances  of  the  basic  chemical  elements: 
carbon,  hydrogen,  oxygen  and  nitrogen  in  Eqs.  (31)— (34): 


<29) 

ln[A  eXP(RS)]  -  +%ne)  =  0  (30) 

fc%-(f°c0hI+fc^+fc0o>)  =°  (31) 

4 FW  +  2 F*V>0  -  (4F°H>  +  2 F°%  +  2F^j))  =  0  (32) 

4®  -  (4 o1  +  Fcol  +  Fh!o)  =  0  (33) 

4®  “  fNa 1  =  0  (34) 


where  F  designates  the  molar  flow  rate  [mol  s_1  ],  superscripts  i  and 
o  represents  input  and  output  flows  and  subscripts  describes  the 
respective  chemical  species.  The  additional  variables  mijne  and  riune 
are  the  Arrhenius  line  coefficients.  All  of  the  equations  presented 
above  are  applied  for  the  calculation  of  each  experimental  point 
(configuration)  -  the  total  number  of  equations  used  in  the 
mathematical  model  equals  6n,  where  n  is  the  amount  of  experi¬ 
mental  points. 

The  total  number  of  measurement  points  used  in  the  analysis 
was  n  =  52.  Table  3  presents  the  vector  containing  the  measured 
and  unknown  variables.  The  beginning  part  contains  measured 
values  (10  variables  for  each  of  the  measurement  points,  in  total 
520  variables)  and  then  the  unknowns.  The  most  important  un¬ 
knowns  are  the  empirical  parameters  that  describe  the  process:  a, 
b,  E  and  A.  Moreover,  the  additional  variables,  like  Arrhenius  line 
coefficients  and  output  flows  of  water,  nitrogen  and  all  compo¬ 
nents,  F^®  ar>d  F^  (defined  for  each  measured  point)  are 
incorporated  into  the  process  model  and  the  total  number  of  un¬ 
knowns  equal  to  162.  Table  3  presents  also  the  assumed  errors  for 
all  the  variables,  which  defines  the  covariance  matrix.  The  un¬ 
certainties  of  the  measured  values  were  determined  on  the  basis  of 
the  error  analysis  conducted  for  the  respective  measuring  devices. 
In  the  case  of  the  unknowns  the  initial  errors  were  assumed  as  100% 
of  the  initiatory  values,  which  were  determined  in  the  preliminary 
calculations. 

After  the  definition  of  all  elements  of  the  mathematical  model 
was  made,  it  was  possible  to  apply  the  GLS  algorithm:  linearize  the 
constraint  equations,  determine  the  Jacobi’s  matrix  and  compute 
the  correction  vectors  VB  and  new  error  bounds  for  the  corrected 
data  (covariance  matrix  Cvb).  The  values  of  the  corrected  unknowns 
and  its  errors  are  presented  in  Table  4,  and  the  corrected  kinetic 
equation  is: 
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in  the  correlation  with  Gas  Hourly  Space  Velocity  (GHSV)  at  the  condition  (SC  =  3  and 
NC  =  3),  (c)  in  the  correlation  with  GHSV  at  the  condition  of  the  reaction  temperature 
at  700[°C], 


Rst  =  Wot-  6.472 

x  10-^exp(=illf^)(pCH,)“»(p„lO)-"“  (35) 

Table  5  presents  the  results  obtained  for  one  of  the  measure¬ 
ment  points,  and  similar  results  can  be  observed  in  the  case  of  all  of 


them.  It  can  be  noted  that  the  corrected  data  has  smaller  uncer¬ 
tainty  than  the  initially  assumed  one.  It  is  an  effect  of  applying 
constraint  equations,  which  provides  additional  information  about 
the  analysed  system.  The  biggest  improvement  can  be  seen  in  the 
case  of  temperature  measurement  as  well  as  in  the  results  obtained 
from  the  gas  chromatography— output  flows  of  CH4,  H2,  CO2  and  CO. 


io-j  +  1 

lOj  +  2 
lOj  +  3 
lOj  +  4 


t « 

3 

Sr 

"C 


[°C] 

[g] 


l-1] 

l-1] 

l-1] 


±0.002-F*®  „ 
±0.04F®Ov 
±0.002.F®„ 


±0.5 

±0.5 


points  -nj  =  0,1... .n-l 


1.554  x  10  3 
117.22  x  103 


6.472  x  10-3 
1213  x  103 


±2.8  x  103 


The  improvement  obtained  by  the  use  of  the  GLS  algorithm  can 
be  seen  clearly  by  comparison  of  the  residuum  of  the  model 
equations  Wj  solved  with  the  initial  and  corrected  data.  For  the 
kinetic  reforming  equation  in  Eq.  (29),  the  sum  of  the  residuum  for 
all  measurement  points  before  the  application  of  the  GLS  algorithm 
is  1.4  x  10-2  and  after  data  correction  it  is  9.798  x  10-4.  Similar 
dependence  can  be  observed  in  the  case  of  other  equations  —  the 
sum  of  the  residuum  for  all  equations  applied  to  all  the  measure¬ 
ment  points  is  11.47  before  the  application  of  the  GLS  algorithm  and 
3.70  after. 


5.  Conclusions 

The  present  paper  highlighted  the  challenge  in  the  study  of 
reaction  kinetics  derived  in  the  methane/steam  reforming  process 
with  a  novel  approach  combining  the  Orthogonal  Least  Squares 
mathematical  methodology.  A  calculation  procedure  incorporating 


Table  5 

Analysed  results  -  variables  for  one  of  the  measurement  point  (initial  and  corrected 


Variables  Initial  value  ±Initial  eri 

-or  Correction 

Corrected 

±Corrected 

T 

700 

±2.49 

-0.004 

700.00 

±0.36 

AP 

FCH4  „ 

Fko, 

5.0  x  10- 
35 

0.13 

3  ±1.50  x  1< 
±0.07 
±5.2  x  10 

r5  7.85  x  10-9 
-4.6  x  10-' 
3  -8.9  x  10"' 

5.00  x  10- 
1  35.000 

1  0.1309 

3  ±1.50  x  10-5 
±0.070 
±5.1  x  10~3 

140 

18.328 

14.553 

±0.28 

±1 

1.11  x  10"1 

0.6171 

-0.4787 

7  140 

17.71 

15.03 

±0.28 

±0.70 

±0.75 

0.532 

±0.5 

-0.1353 

0.67 

±0.38 

the  Generalized  Least  Squares  (GLS)  method  has  been  proposed  to 
estimate  the  basic  parameters  of  the  methane/steam  reforming 
kinetics  over  the  Ni/YSZ,  typically  used  for  SOFC  anodic  catalytic 
material.  The  introduced  mathematical  model  has  been  based  only 
on  the  fundamental  physical  equations  describing  the  chemical 
process,  and  it  can  be  noted  that  it  does  not  contain  any  additional 
assumptions  in  accordance  with  the  precise  stoichiometry  of  the 
occurring  reactions.  This  property  provides  the  greater  versatility 
to  the  proposed  methodology  in  comparison  with  approaches  used 
in  the  previous  studies. 

Three  main  phases  of  the  conducted  research  can  be  pointed: 
the  experimental  investigation  of  the  methane/steam  reforming 
process  over  a  fine  powder  catalyst  of  NiO/YSZ  after  reduction,  the 
numerical  analysis  of  collected  data  which  aimed  to  find  the  initial 
approximations  of  empirical  parameters  describing  the  process, 
and  finally  the  application  of  the  GLS  algorithm. 

The  method  proposed  in  the  present  paper  not  only  allows  to 
obtain  the  most  probable  values  of  the  estimated  parameters  but 
also  makes  possible  the  a  posteriori  evaluation  of  the  errors  of 
directly  measured  variables  and  unknown  parameters.  It  was 
demonstrated  that  the  GLS  method  is  useful  in  securing  higher 
accuracy  of  measured  data  and  decreasing  the  value  of  residuum  of 
constraint  equations.  The  present  results  primarily  proved  an 
applicability  of  the  GLS  methodology  into  the  investigation  of  a 
chemical  reaction  process  for  deriving  its  reaction  kinetics.  More¬ 
over,  it  can  be  noted  that  the  GLS  method  has  the  potential  to 
provide  the  objective  criteria  for  the  formal  evaluation  and  falsifi¬ 
cation  of  different  mathematical  models  of  the  methane  reforming 
process.  Furthermore,  the  proposed  algorithm  can  be  easily 
adapted  to  the  mathematical  modelling  of  the  fuel  reforming  pro¬ 
cess  for  different  catalysts  and  experimental  conditions.  This 
approach  can  clarify  the  divergence  between  the  published  con¬ 
tradictory  results  about  the  reaction  kinetics’  parameters.  The 
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authors  of  the  present  paper  strongly  believe  that  the  application  of 
the  GLS  method  can  provide  a  great  benchmark  in  the  evaluation  of 
the  reaction  process  and  qualify  the  level  of  mathematical  model¬ 
ling  of  reaction  processes. 
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